%% Assignment 8.3
% CIE4440

close all;clearvars;clc;

%% initiate data
plot    = 0;

dx      = 1;
x       = 0:dx:3000;
width_he= 8;
startDepth=2.35;
g       = 9.81;

QSmall  = 1;
QLarge  = 40;
QJoint  = QSmall+QLarge;

ChSmall = 20;
ChLarge = 20;
ChJoint = 20;

ibSmall = 1/1000;
ibLarge = 1/1000;
ibJoint = 1/1000;

%% 8.3.1 How can for the assumption B/ he = 8 an equilibrium depth be
% calculated from discharge Q, roughness coefficient C, and slope i.

% This is a rewrite of the chezy equation for equilibrium depth.
heSmall = (QSmall^2/(ChSmall^2*width_he^2^2*ibSmall))^.2;
heLarge = (QLarge^2/(ChLarge^2*width_he^2^2*ibLarge))^.2;
heJoint = (QJoint^2/(ChJoint^2*width_he^2^2*ibJoint))^.2;



%% 8.3.2 We wish to investigate how far upstream of the confluence a
% gauging station has to be constructed in order to have no influence from
% backwater. Assume an accuracy of 1% as acceptable. Investigate the required
% distance upstream from the confluence for variable discharges in each of
% the branches ranching from 40 m3/sec to 5 m3/sec. with steps of 5 m3/sec
% and between 5 m3/sec and 1 m3/sec with steps of 1 m3/sec . Present the
% results in a table for each combination of discharges.

if plot
    figure(3)
    cax     = gca;
    set(cax,'NextPlot','add');
    title('BWC');
    xlabel('x m');
    ylabel('h m');
end

[QSmall,QLarge]=meshgrid(unique([1:5 5:5:40]));

for ii=1:numel(QSmall)
    QJoint(ii)  = QSmall(ii)+QLarge(ii);
    
    heSmall(ii) = (QSmall(ii)^2/(ChSmall^2*width_he^2*ibSmall))^0.2;
    heLarge(ii) = (QLarge(ii)^2/(ChLarge^2*width_he^2*ibLarge))^0.2;
    heJoint(ii) = (QJoint(ii)^2/(ChJoint^2*width_he^2*ibJoint))^0.2;
    
    % Where does the to the power 1/3 comes from? Fr=1
    hcSmall(ii) = (QSmall(ii)^2/(g*width_he^2))^(1/3);
    hcLarge(ii) = (QLarge(ii)^2/(g*width_he^2))^(1/3);
    
    for i=1:length(x)
        if i==1
            hSmall(i)=heJoint(ii);
            hLarge(i)=heJoint(ii);
        else
            hSmall(i)=hSmall(i-1)+dhdxSmall(i-1)*dx;
            hLarge(i)=hLarge(i-1)+dhdxLarge(i-1)*dx;
        end
        
        dhdxSmall(i)=-ibSmall*(hSmall(i)^3-heSmall(ii)^3)/(hSmall(i)^3-hcSmall(ii)^3);
        dhdxLarge(i)=-ibLarge*(hLarge(i)^3-heLarge(ii)^3)/(hLarge(i)^3-hcLarge(ii)^3);
        
        h_heSmall(i)=(hSmall(i)/heSmall(ii))-1;
        h_heLarge(i)=(hLarge(i)/heLarge(ii))-1;
        
    end
    %     fprintf('%f,%f\n',QSmall(ii),QLarge(ii));
    
    placeGaugeSmall(ii)=x(find(h_heSmall<0.01,1,'first'));
    placeGaugeLarge(ii)=x(find(h_heLarge<0.01,1,'first'));
    
    if plot
        plot(x,hSmall,x,hLarge)
    end
    
end
tableArray=[QSmall(:) QLarge(:) placeGaugeSmall(:) placeGaugeLarge(:)];

fprintf('Q Small | Q Large | Location gauge small | Location gauge large\n')
disp(tableArray);